Quais variáveis escolher para o meu modelo?
Edson Silva-Júnior
Uma análise sobre os principais métodos de escolha de variáveis de modelos, com aplicações em R.
A Ecologia é a ciência responsável por tentar entender como os organismos reagem com a interação com o seu ambiente, incluindo outros organismos. Por causa disso, a ecologia utiliza de forte uso de modelagens matemáticas e estatísticas para medir se as relações pesquisadas existem e o quanto elas são fortes. Então como a ecologia trabalha em modelar relações dos seres vivos com o mundo, é importante que seus modelos estejam bem ajustados, minimizando os erros associados. Uma das etapas mais importantes da criação de um modelo é a escolha das principais variáveis utilizadas no mesmo. Aí entramos em uma das maiores perguntas quando se fala de modelos: quais são as variáveis que eu devo utilizar no meu modelo?
O que são modelos?
A terceira década do século XXI (2021-2030) será marcada pelo uso de tecnologias que “pensam por si só”. Nesse época se fala muito dos modelos generativos de informação, como o famoso ChatGPT. E por mais que engenheiros de softwares voltados para Inteligências Artificiais e analistas de dados econômicos, com uma certa frequêncai, tentem reinvidicar que seus modelos são perfeitos, a realidade é que todos estão errados. Mas afinal, o que são modelos?
Modelos, no contexto matemático e estatístico, são
tentativas de representar fenômenos do mundo real para números
e algorítimos, que gerem resultados que possam ser interpretados como a
expressão destes fenômenos. Mas a palavra chave é exatamente essa:
tentativas. Acontece que nenhum modelo consegue ser o fenômeno
modelado de fato, mas sim partes deles. Alguns mais parecidos, outros
menos. Mas, no fim, nenhum representa de fato os fenômenos estudados. A
diferença entre os resultados que um modelo encontrou para o fenômeno do
mundo real que este tentou modelar é o que chamamos de
erro. Logo, quanto menos distância existe entre os
resultados de um modelo e o que aquele fenômeno no mundo real é, menos
errado o modelo em questão é. Dentro da ecologia, um exemplo clássico
disso são os modelos de crescimento populacional, que não são perfeitas,
possuem erros, mas são extremamente úteis nas tomadas de decisão
envolvendo conservação de biodiversidade e controle de epidemias.
Existem diferentes formas de avaliar os erros de um modelo. Algumas das principais formas empregadas na Ecologia são:
A forma mais simples, que é apenas olhar para o que o modelo gerou e comparar o quão diferente isto é do fenômeno estudado;
Com estatística descritiva, onde os erros podem ser medidos através de avaliações com medidas de tendência central, como a média (\(\bar{x}\), \(\mu\)), aliadas à medidas de dispersão, como o desvio padrão (\(s\), \(\sigma\));
Com estatística frequentista, através de testes de hipótese, os erros podem ser avaiados como a diferença entre a estatística calculada (testes T, Z, W, F, X2, etc.) e a estatística crítica;
Também relacionada à estatística frequentista, através de coeficientes, como o coeficiente de correlação (r) e o coeficiente de determinação (R2), e resultados sobre significância daqueles resultados ao acaso (valor de p);
Através de redução de dimensionalidade das variáveis, como na estatística multidimensional, a Análise de Componentes Principais (Principal Components Analysis - PCA), e a comparação dos autovalores das dimensões reduzidas;
Qualquer outra forma matemática, estatística e/ou teórico-conceitual que permita medir o grau de falha de um modelo em medir um fenômeno de forma perfeita.
É importante destacar que a única coisa que reúne todas essas coisas,
de diferentes tipos de modelos, interpretados como erros é
o fato de que são os resultados de falhas na explição
perfeita de um fenômeno do mundo real. Mas, de forma geral,
diferentes tipos de modelos apresentam diferentes formas de expressar
seus erros, então a forma de avaliar os erros de um modelo precisa ser
adequada não apenas para as naturezas matemáticas do modelo avaliados,
mas também para a teoria biológica por trás da medição daquele modelo.
Um mesmo modelo pode estar mais errado ou menos errado, com os mesmos
dados e sem alterar seus parâmetros, a depender da forma e das
ferramentas com o qual este foi avaliado.
E no que exatamente isso afeta na escolha de variáveis para um
modelo? Adicionalmente ao que foi citado no tópico 1 da
lista anterior, a escolha de variáveis pouco relevantes para um modelo
pode fazer com que seu erro aumente, já que estas não apresentam de fato
importância real para o fenômeno estudado, por diferentes motivos. Desta
forma, a escolha das melhores variáveis de um modelo também influencia
no seu desempenho.
A seguir, veremos as principais formas utilizadas dentro da ecologia para a escolha de variáveis, abordando a teoria e a prática utilizando linguagem de programação R.
Importância ecológica e biológica
Antes de qualquer análise estatística, é importante que suas
variáveis façam sentido ecológico e biológico.
Uma confusão relativamente frequente quando se trabalha com ciência de dados é coletar o máximo de váriáveis possíveis nos modelos e assumir que as que apresentaram os maiores desempenhos são os fatores com maior influências para o fenômeno estudado. O grande erro disso é que qualquer uma dessas análises só medem relações matemáticas, mas isso não significa que essas relaçõees de fato existem. Para exemplificar de forma mais clara esses efeitos, um site ótimo é o Spurious Correlations (Correlações Espúrias, em inglês): esse site mostra como coisas sem correlação alguma no mundo real podem apresentar correlações matemáticas dos seus dados, ou, como seu nome sugere, correlações espúrias.
Esses problemas se tornam mais graves quando lembramos que, em contextos ecológicos, variáveis ambientais (sejam bióticas ou abióticas) podem ou não ter influência, dependendo do objeto de estudo. A título de exemplo, Russildi et al. (2016), ao avaliarem os fatores que mais influenciavam a diversidade de anfíbios e répteis em florestas tropicais no México, encontraram que a diversidade de anfíbios e répteis são ambos influenciados por maiores precipitações, mas anfíbios sofreram mais influência de temperaturas menores e maior quantidade de corpos d’água próximos, enquanto que répteis sofreram maior influência das formas dos fragmentos e maiores temperaturas. Note que é o mesmo recorte geográfico e temporal, mas diferentes grupos sofreram influências de diferentes variáveis.
Por causa disso, o primeiro passo na construção de um modelo
útil é conhecer a ecologia e a biologia do objeto estudado.
Conhecer a biologia, a história natural, a ecofisiologia, interações
bióticas e quaisquer informações ecológicas relevantes do seu objeto de
estudo é fundamental para que seu modelo esteja menos errado,
minimizando esforços e recursos desperdiçados com variáveis pouco
interessantes para responder a pergunta de pesquisa em questão.
Multicolinearidade
A Multicolinearidade pode ser definida quando diferentes variáveis
variam de foram semelhante, ou sejam:
apresentam linearidades semelhantes. Nos estudos
ecológicos, isso apresenta dois grandes problemas:
Na natureza, é extremamente comum que fatores ambientais apresentem correção, e então sofram de multicolinearidade. Um exemplo clássico disso é a temperatura e a umidade de um local: na maioria dos ecossistemas terrestres, sobretudo nos trópicos, é normal que áreas com maiores taxas de umidade apresentem menor temperatura, e vice-versa. Outro exemplo comum desses fenômenos em ecossistemas terrestres tropicais é a correlação entre diversidade de organismos saxicólas (organismos que vivem em serrapilheira) e o grau de abertura de dossel;
Diversos tipos de modelos utilizados para analisar fenômenos ecológicos partem da premissa de independência das unidades amostrais: cada ponto, local e repetição dos dados precisam ser informações novas, que os seus valores não dependam dos outros valores dos outros pontos, mesmo que estes apresentem valores semelhantes. Quando essa premissa não é obedecida, o modelo se torna inflado, e então este acaba gerando mais influência para uma tendência de dados não por ela, biologica, ecológica e matematicamente, fazer mais sentido, mas sim pelo fato desta informação se tornar repetida.
Devido a estes dois principais fatores, um parte importantes dos modelos ecológicos envolvem não apenas detectar, mas lidar com a multicolinearidade das variáveis.
Remoção de multicolinearidade
A técnica mais simples envolvendo lidar com multicolinearidade é a remoção de variáveis que apresentam correlação, selecionando, de forma manual e baseado na biologia e ecologia do objeto de estudo, as variáveis mais interessantes para responder a pergunta ecológica.
Matriz de correlação
Uma forma simples de detectar e medir a correlação entre variáveis de
um modelo é através dos índices de correlação: índices
matemáticos que calculam não apenas se duas ou mais variáveis variam de
forma semelhante, mas o quanto essa correlção é forte.
Existem diferentes métodos para isso, mas os dois principais dentro
de estudos ecológicos são os Índice de Pearson e
Índice de Spearman: o Índice de Pearson mede correlações
lineares, enquanto que o índice de Spearman mede correlações não
lineares. Estes índices variam de -1 a 1, funcionando da seguinte
forma:
r = 1: correlação diretamente proporcional perfeita.
r = 0: ausência de correlação;
r = -1: correlação inversamente proporcional perfeita.
A escolhas entre os índices de Pearson e Spearman é uma discussão bastante antiga dentro da estatística, com diversos campos de estudo apontando motivos para escolher entre um e outro. Dentro da ecologia, quando se deseja identificar e medir multicolinearidade entre variáveis ambientais, a maioria dos pesquisadores empregam o Índice de Spearman, divido a relação entre diversas variáveis ambientais não ocorrer de forma linear. Contudo, é importante que a escolha entre esses índices sejam pensado para cada estudo como um caso em particular
Para o nosso exemplo, imagine a seguinte situação: você está estudando a relação entre a taxa de herbivoria que três espécies de flores do gênero Iris (I. setosa, I. versicolor, I. virginica) sofrem e a morfometria de suas pétalas e sápalas, tendo dados do comprimento e largura das pétalas e sépalas dos indivíduos das três espécies. Antes de criar um modelo você irá testar a multicolinearidade através de índices de correlação. Para esse exemplo, serão usados tanto o índice de Pearson quando de Spearman.
Para o código em R, usaremos o dataset iris, nativo do
R, e os pacotes tidyverse, para o tratamento e a
visualização de dados e reshape2, para a transformação da
matriz de correlação.
Primeiro, vamos carregar nossos pacotes. Certifique-se que estão instalados antes de começar.
# Carregando os pacotes
library(tidyverse)## ── Attaching core tidyverse packages ────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ──────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(reshape2)##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smithsEm seguida, vamos carregar nossa base de dados, através da função
data(), informando o nome da base de dados entre aspas, e
então chamando-a.
# Carreando a base de dados
data("iris")
# Visualizando a base de dados
iris## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## 11 5.4 3.7 1.5 0.2 setosa
## 12 4.8 3.4 1.6 0.2 setosa
## 13 4.8 3.0 1.4 0.1 setosa
## 14 4.3 3.0 1.1 0.1 setosa
## 15 5.8 4.0 1.2 0.2 setosa
## 16 5.7 4.4 1.5 0.4 setosa
## 17 5.4 3.9 1.3 0.4 setosa
## 18 5.1 3.5 1.4 0.3 setosa
## 19 5.7 3.8 1.7 0.3 setosa
## 20 5.1 3.8 1.5 0.3 setosa
## 21 5.4 3.4 1.7 0.2 setosa
## 22 5.1 3.7 1.5 0.4 setosa
## 23 4.6 3.6 1.0 0.2 setosa
## 24 5.1 3.3 1.7 0.5 setosa
## 25 4.8 3.4 1.9 0.2 setosa
## 26 5.0 3.0 1.6 0.2 setosa
## 27 5.0 3.4 1.6 0.4 setosa
## 28 5.2 3.5 1.5 0.2 setosa
## 29 5.2 3.4 1.4 0.2 setosa
## 30 4.7 3.2 1.6 0.2 setosa
## 31 4.8 3.1 1.6 0.2 setosa
## 32 5.4 3.4 1.5 0.4 setosa
## 33 5.2 4.1 1.5 0.1 setosa
## 34 5.5 4.2 1.4 0.2 setosa
## 35 4.9 3.1 1.5 0.2 setosa
## 36 5.0 3.2 1.2 0.2 setosa
## 37 5.5 3.5 1.3 0.2 setosa
## 38 4.9 3.6 1.4 0.1 setosa
## 39 4.4 3.0 1.3 0.2 setosa
## 40 5.1 3.4 1.5 0.2 setosa
## 41 5.0 3.5 1.3 0.3 setosa
## 42 4.5 2.3 1.3 0.3 setosa
## 43 4.4 3.2 1.3 0.2 setosa
## 44 5.0 3.5 1.6 0.6 setosa
## 45 5.1 3.8 1.9 0.4 setosa
## 46 4.8 3.0 1.4 0.3 setosa
## 47 5.1 3.8 1.6 0.2 setosa
## 48 4.6 3.2 1.4 0.2 setosa
## 49 5.3 3.7 1.5 0.2 setosa
## 50 5.0 3.3 1.4 0.2 setosa
## 51 7.0 3.2 4.7 1.4 versicolor
## 52 6.4 3.2 4.5 1.5 versicolor
## 53 6.9 3.1 4.9 1.5 versicolor
## 54 5.5 2.3 4.0 1.3 versicolor
## 55 6.5 2.8 4.6 1.5 versicolor
## 56 5.7 2.8 4.5 1.3 versicolor
## 57 6.3 3.3 4.7 1.6 versicolor
## 58 4.9 2.4 3.3 1.0 versicolor
## 59 6.6 2.9 4.6 1.3 versicolor
## 60 5.2 2.7 3.9 1.4 versicolor
## 61 5.0 2.0 3.5 1.0 versicolor
## 62 5.9 3.0 4.2 1.5 versicolor
## 63 6.0 2.2 4.0 1.0 versicolor
## 64 6.1 2.9 4.7 1.4 versicolor
## 65 5.6 2.9 3.6 1.3 versicolor
## 66 6.7 3.1 4.4 1.4 versicolor
## 67 5.6 3.0 4.5 1.5 versicolor
## 68 5.8 2.7 4.1 1.0 versicolor
## 69 6.2 2.2 4.5 1.5 versicolor
## 70 5.6 2.5 3.9 1.1 versicolor
## 71 5.9 3.2 4.8 1.8 versicolor
## 72 6.1 2.8 4.0 1.3 versicolor
## 73 6.3 2.5 4.9 1.5 versicolor
## 74 6.1 2.8 4.7 1.2 versicolor
## 75 6.4 2.9 4.3 1.3 versicolor
## 76 6.6 3.0 4.4 1.4 versicolor
## 77 6.8 2.8 4.8 1.4 versicolor
## 78 6.7 3.0 5.0 1.7 versicolor
## 79 6.0 2.9 4.5 1.5 versicolor
## 80 5.7 2.6 3.5 1.0 versicolor
## 81 5.5 2.4 3.8 1.1 versicolor
## 82 5.5 2.4 3.7 1.0 versicolor
## 83 5.8 2.7 3.9 1.2 versicolor
## 84 6.0 2.7 5.1 1.6 versicolor
## 85 5.4 3.0 4.5 1.5 versicolor
## 86 6.0 3.4 4.5 1.6 versicolor
## 87 6.7 3.1 4.7 1.5 versicolor
## 88 6.3 2.3 4.4 1.3 versicolor
## 89 5.6 3.0 4.1 1.3 versicolor
## 90 5.5 2.5 4.0 1.3 versicolor
## 91 5.5 2.6 4.4 1.2 versicolor
## 92 6.1 3.0 4.6 1.4 versicolor
## 93 5.8 2.6 4.0 1.2 versicolor
## 94 5.0 2.3 3.3 1.0 versicolor
## 95 5.6 2.7 4.2 1.3 versicolor
## 96 5.7 3.0 4.2 1.2 versicolor
## 97 5.7 2.9 4.2 1.3 versicolor
## 98 6.2 2.9 4.3 1.3 versicolor
## 99 5.1 2.5 3.0 1.1 versicolor
## 100 5.7 2.8 4.1 1.3 versicolor
## 101 6.3 3.3 6.0 2.5 virginica
## 102 5.8 2.7 5.1 1.9 virginica
## 103 7.1 3.0 5.9 2.1 virginica
## 104 6.3 2.9 5.6 1.8 virginica
## 105 6.5 3.0 5.8 2.2 virginica
## 106 7.6 3.0 6.6 2.1 virginica
## 107 4.9 2.5 4.5 1.7 virginica
## 108 7.3 2.9 6.3 1.8 virginica
## 109 6.7 2.5 5.8 1.8 virginica
## 110 7.2 3.6 6.1 2.5 virginica
## 111 6.5 3.2 5.1 2.0 virginica
## 112 6.4 2.7 5.3 1.9 virginica
## 113 6.8 3.0 5.5 2.1 virginica
## 114 5.7 2.5 5.0 2.0 virginica
## 115 5.8 2.8 5.1 2.4 virginica
## 116 6.4 3.2 5.3 2.3 virginica
## 117 6.5 3.0 5.5 1.8 virginica
## 118 7.7 3.8 6.7 2.2 virginica
## 119 7.7 2.6 6.9 2.3 virginica
## 120 6.0 2.2 5.0 1.5 virginica
## 121 6.9 3.2 5.7 2.3 virginica
## 122 5.6 2.8 4.9 2.0 virginica
## 123 7.7 2.8 6.7 2.0 virginica
## 124 6.3 2.7 4.9 1.8 virginica
## 125 6.7 3.3 5.7 2.1 virginica
## 126 7.2 3.2 6.0 1.8 virginica
## 127 6.2 2.8 4.8 1.8 virginica
## 128 6.1 3.0 4.9 1.8 virginica
## 129 6.4 2.8 5.6 2.1 virginica
## 130 7.2 3.0 5.8 1.6 virginica
## 131 7.4 2.8 6.1 1.9 virginica
## 132 7.9 3.8 6.4 2.0 virginica
## 133 6.4 2.8 5.6 2.2 virginica
## 134 6.3 2.8 5.1 1.5 virginica
## 135 6.1 2.6 5.6 1.4 virginica
## 136 7.7 3.0 6.1 2.3 virginica
## 137 6.3 3.4 5.6 2.4 virginica
## 138 6.4 3.1 5.5 1.8 virginica
## 139 6.0 3.0 4.8 1.8 virginica
## 140 6.9 3.1 5.4 2.1 virginica
## 141 6.7 3.1 5.6 2.4 virginica
## 142 6.9 3.1 5.1 2.3 virginica
## 143 5.8 2.7 5.1 1.9 virginica
## 144 6.8 3.2 5.9 2.3 virginica
## 145 6.7 3.3 5.7 2.5 virginica
## 146 6.7 3.0 5.2 2.3 virginica
## 147 6.3 2.5 5.0 1.9 virginica
## 148 6.5 3.0 5.2 2.0 virginica
## 149 6.2 3.4 5.4 2.3 virginica
## 150 5.9 3.0 5.1 1.8 virginicaPara calcularmos a multicolinearidade, primeiro temos que selecionar
apenas as variáveis numéricas, com as funções
dplyr::select(), para selecionar colunas, e
dplyr::where(), para informar que apenas serão selecionadas
as variáveis numéricas (is.numeric), ambas do pacote
dplyr, que faz parte dos pacotes carregados pelo
tidyverse. Em seguida, utilizaremos a função nativa
cor(), informando que será utilizado o índice de Pearson
(method = "pearson"). Tudo isso atribuido (operador
<-) ao objeto cor_pearson
Note que as funções estão se operando uma após a outra, ao invés de
todas estarem aninhadas umas dentros das outras. Isso significa que uma
função está recendo os dados gerados pela função anterior. Isso é
possível através do pipe (%>%), que é ativado quando o
pacote tidyverse é ativado. Ele permite deixar o código
mais legível e menos confuso.
# Correlação de Pearson
cor_pearson <- iris %>%
dplyr::select(dplyr::where(is.numeric)) %>% # selecionando apenas as variáveis com dados numéricos
cor(method = "pearson") # especificando que o índice é o de Pearson
## Matriz de correlação
cor_pearson## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
## Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
## Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
## Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000Note que a variável do comprimento da sépala
(Sepal.Length) está altamente correlacionada com as
variáveis de comprimento e largura das pétalas
(Petal.Length e Petal.Width, respectivamente),
devido seu alto valor de correlação (r > 0.7). Para ajudar na
interpretação, podemos criar um gráfico para visualizar estas relações,
através de ggplot.
## Gráfico
cor_pearson[upper.tri(cor_pearson)] <- NA # removendo a matriz triangular superior através de transforma-la em NA
cor_pearson %>%
reshape2::melt() %>% # ajustando a matriz para um formato de data frame
na.omit() %>% # removendo os NAs gerados
ggplot(aes(Var1, Var2, fill = value)) +
geom_tile(color = "black") +
geom_text(aes(label = value %>% round(2))) +
labs(fill = "Índice de Correllação") +
scale_fill_viridis_c(limits = c(-1, 1), option = "turbo", breaks = seq(-1, 1, 0.2)) +
guides(fill = guide_colorbar(title.position = "top",
title.hjust = 0.5,
barwidth = 25,
barheight = 1.5,
frame.colour = "black",
ticks.colour = "black")) +
theme(axis.text = ggplot2::element_text(color = "black", size = 10),
axis.text.x = ggplot2::element_text(color = "black", size = 10, angle = 90, vjust = 1),
axis.title = element_blank(),
strip.background = ggplot2::element_rect(color = "black"),
axis.line = ggplot2::element_line(color = "black"),
panel.background = ggplot2::element_rect(color = "black", fill = "white"),
legend.position = "bottom",
legend.title = element_text(color = "black", size = 12))
Vamos então repetir o mesmo processo, mas agora para Spearman.
# Correlação de Spearman
cor_spearman <- iris |>
dplyr::select(dplyr::where(is.numeric)) |> # selecionando apenas as variáveis com dados numéricos
cor(method = "spearman")
## Matriz de correlação
cor_spearman## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000000 -0.1667777 0.8818981 0.8342888
## Sepal.Width -0.1667777 1.0000000 -0.3096351 -0.2890317
## Petal.Length 0.8818981 -0.3096351 1.0000000 0.9376668
## Petal.Width 0.8342888 -0.2890317 0.9376668 1.0000000## Gráfico
cor_spearman[upper.tri(cor_spearman)] <- NA
cor_spearman %>%
reshape2::melt() %>%
na.omit() %>%
ggplot(aes(Var1, Var2, fill = value)) +
geom_tile(color = "black") +
geom_text(aes(label = value %>% round(2))) +
labs(fill = "Índice de Correllação") +
scale_fill_viridis_c(limits = c(-1, 1), option = "turbo", breaks = seq(-1, 1, 0.2)) +
guides(fill = guide_colorbar(title.position = "top",
title.hjust = 0.5,
barwidth = 25,
barheight = 1.5,
frame.colour = "black",
ticks.colour = "black")) +
theme(axis.text = ggplot2::element_text(color = "black", size = 10),
axis.text.x = ggplot2::element_text(color = "black", size = 10, angle = 90, vjust = 1),
axis.title = element_blank(),
strip.background = ggplot2::element_rect(color = "black"),
axis.line = ggplot2::element_line(color = "black"),
panel.background = ggplot2::element_rect(color = "black", fill = "white"),
legend.position = "bottom",
legend.title = element_text(color = "black", size = 12))
Note que, tanto para Pearson quanto Spearman, a combinação sem multicolinearidade é o comprimento e largura da sépala, se seguirmos o nosso ponto de corte como r > |0.7|.
Na matemática, usamos a expressão | | para informar que
o número pode assumir tantos valores positivos quanto negativos. Ou seja
r > |0.7| = r > 0.7 ou r < -0.7.
Também vale destacar que esse valor é arbitário, e qualquer outro valor, desde que faça sentido, poderia ser utilizado.
Fatores de Inflação da Variância (Variance Inflaction Factors - VIF)
Outra forma de detectar a multicolinearidade é através da criação prévia dos modelos, e então os ajustar. Isso pode ser uma vantagem, principalmente, porque, mesmo após a remoção de multicolinearidade através da remoção de variáveis correlacionadas, é possível que o modelo esteja inflado, devido à natureza matemática das variáveis. Por causa disso, uma alternativa é análise de Fatores de inflação da Variância (Variance Inflaction Factors - VIF).
o VIF é uma análise que avalia quantififando o quanto o erro padrão dos coeficientes estimados (\(\beta\) de uma regressão, por exemplo) se tornaram inflados, devido à multicolinearidade das variáveis do modelo. A vantagem deste método é que ele avalia o modelo pronto, e não apenas as suas variáveis, informando se o modelo precisa ou não que alguma variável seja removida.
Em modelos que testam se uma ou mais variáveis influenciam o comportamento de um grupo de dados, como na regressão linear, coeficientes estimados são uma estimativa média do quanto a variação dos preditores de um modelo afetam a variação dos dados das variáveis respostas. Por exemplo, se criarmos um modelo para estimar o quanto a predação (incidência de predadores) afeta a frequência de cantos de anfíbios (Hz), e o coeficiente estimado for -15,5, isso significa que para cada aumento em 1 nos valores dos dados de predação, os valores dos dados de frequência de canto diminuem em 15,5.
Contudo, essa relação não ocorre de forma perfeita. No nosso exemplo,
diferentes indivíduos de anfíbios, mesmo da mesma espécie, responderão
de forma de diferente à predação, devido a fatores como tamanho do
corpo, morfomoetria de membros e orgãos, ecofisioloigia, estado de
nutrição e hidratação, genética e até personalidade individual. Por isso
os coeficientes estimados são médias, não valores
absolutos. Assim, algumas unidades amostrais serão mais afetadas e
outras menos afetadas, alguns sendo afetadas para diminuir -16,7 e
outras a diminuir -11,63. Por isso, é utilizado a estimativa do erro
padrão, em conjunto com os coeficientes estimados (\(\mu \pm S \bar{x}\)), já que um erro padrão
baixo significa que que os valores do quanto as unidades amostrais são
afetadas orbitam muito próximos da média, o que significa, no nosso
exemplo, que a relação de influência da predação sobre a frequência do
canto dos anfíbios é constante e sofre poucas variância. No nosso
exemplo, um coeficiente de determinação com erro padrão = \(15,5 \pm 0,89\) representa uma grande
constância da influência da predação sobre a redução das frequências de
canto dos anfíbios.
Todas as variáveis de um modelo possuem um valor de VIF associado, o que significa que as variáveis que serão removidas não são aquelas que apresentam um valor discriminante, mas sim que estão fora de um intervalo considerado “permitido”. Alguns autores sugerem que o limite seja de 10, então variáveis com VIF > 10 seriam removidas, enquanto outros autores, mais conservadores, recomendam que seja valores como 5, outro 3 e alguns até 2. Em estudos ecológicos, é comum que o valores de cortes estejam entre 3 e 4 (Zuur et al 2010, Dormann et al 2013)
Vamos continuar o nosso exemplo das flores do gênero Iris,
mas agora vamos simular dados de herbivoria, e criar um modelo linear.
Para isso, vamos criar mais uma variável da taxa de herbivoria, através
da função dplyr::mutate(), do pacote dplyr, e
da função sample(), que é nativa do R, atribuindo isso ao
objeto iris_pred.
set.seed(3214) # controlando a leatoriedade, assim sempre será a mesma sequência aleatória de dados
iris_pred <- iris %>%
dplyr::mutate(Predação = sample(seq(0, 1, length.out = 10000000) %>% round(3), size = 150, replace = TRUE))
iris_pred## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## 11 5.4 3.7 1.5 0.2 setosa
## 12 4.8 3.4 1.6 0.2 setosa
## 13 4.8 3.0 1.4 0.1 setosa
## 14 4.3 3.0 1.1 0.1 setosa
## 15 5.8 4.0 1.2 0.2 setosa
## 16 5.7 4.4 1.5 0.4 setosa
## 17 5.4 3.9 1.3 0.4 setosa
## 18 5.1 3.5 1.4 0.3 setosa
## 19 5.7 3.8 1.7 0.3 setosa
## 20 5.1 3.8 1.5 0.3 setosa
## 21 5.4 3.4 1.7 0.2 setosa
## 22 5.1 3.7 1.5 0.4 setosa
## 23 4.6 3.6 1.0 0.2 setosa
## 24 5.1 3.3 1.7 0.5 setosa
## 25 4.8 3.4 1.9 0.2 setosa
## 26 5.0 3.0 1.6 0.2 setosa
## 27 5.0 3.4 1.6 0.4 setosa
## 28 5.2 3.5 1.5 0.2 setosa
## 29 5.2 3.4 1.4 0.2 setosa
## 30 4.7 3.2 1.6 0.2 setosa
## 31 4.8 3.1 1.6 0.2 setosa
## 32 5.4 3.4 1.5 0.4 setosa
## 33 5.2 4.1 1.5 0.1 setosa
## 34 5.5 4.2 1.4 0.2 setosa
## 35 4.9 3.1 1.5 0.2 setosa
## 36 5.0 3.2 1.2 0.2 setosa
## 37 5.5 3.5 1.3 0.2 setosa
## 38 4.9 3.6 1.4 0.1 setosa
## 39 4.4 3.0 1.3 0.2 setosa
## 40 5.1 3.4 1.5 0.2 setosa
## 41 5.0 3.5 1.3 0.3 setosa
## 42 4.5 2.3 1.3 0.3 setosa
## 43 4.4 3.2 1.3 0.2 setosa
## 44 5.0 3.5 1.6 0.6 setosa
## 45 5.1 3.8 1.9 0.4 setosa
## 46 4.8 3.0 1.4 0.3 setosa
## 47 5.1 3.8 1.6 0.2 setosa
## 48 4.6 3.2 1.4 0.2 setosa
## 49 5.3 3.7 1.5 0.2 setosa
## 50 5.0 3.3 1.4 0.2 setosa
## 51 7.0 3.2 4.7 1.4 versicolor
## 52 6.4 3.2 4.5 1.5 versicolor
## 53 6.9 3.1 4.9 1.5 versicolor
## 54 5.5 2.3 4.0 1.3 versicolor
## 55 6.5 2.8 4.6 1.5 versicolor
## 56 5.7 2.8 4.5 1.3 versicolor
## 57 6.3 3.3 4.7 1.6 versicolor
## 58 4.9 2.4 3.3 1.0 versicolor
## 59 6.6 2.9 4.6 1.3 versicolor
## 60 5.2 2.7 3.9 1.4 versicolor
## 61 5.0 2.0 3.5 1.0 versicolor
## 62 5.9 3.0 4.2 1.5 versicolor
## 63 6.0 2.2 4.0 1.0 versicolor
## 64 6.1 2.9 4.7 1.4 versicolor
## 65 5.6 2.9 3.6 1.3 versicolor
## 66 6.7 3.1 4.4 1.4 versicolor
## 67 5.6 3.0 4.5 1.5 versicolor
## 68 5.8 2.7 4.1 1.0 versicolor
## 69 6.2 2.2 4.5 1.5 versicolor
## 70 5.6 2.5 3.9 1.1 versicolor
## 71 5.9 3.2 4.8 1.8 versicolor
## 72 6.1 2.8 4.0 1.3 versicolor
## 73 6.3 2.5 4.9 1.5 versicolor
## 74 6.1 2.8 4.7 1.2 versicolor
## 75 6.4 2.9 4.3 1.3 versicolor
## 76 6.6 3.0 4.4 1.4 versicolor
## 77 6.8 2.8 4.8 1.4 versicolor
## 78 6.7 3.0 5.0 1.7 versicolor
## 79 6.0 2.9 4.5 1.5 versicolor
## 80 5.7 2.6 3.5 1.0 versicolor
## 81 5.5 2.4 3.8 1.1 versicolor
## 82 5.5 2.4 3.7 1.0 versicolor
## 83 5.8 2.7 3.9 1.2 versicolor
## 84 6.0 2.7 5.1 1.6 versicolor
## 85 5.4 3.0 4.5 1.5 versicolor
## 86 6.0 3.4 4.5 1.6 versicolor
## 87 6.7 3.1 4.7 1.5 versicolor
## 88 6.3 2.3 4.4 1.3 versicolor
## 89 5.6 3.0 4.1 1.3 versicolor
## 90 5.5 2.5 4.0 1.3 versicolor
## 91 5.5 2.6 4.4 1.2 versicolor
## 92 6.1 3.0 4.6 1.4 versicolor
## 93 5.8 2.6 4.0 1.2 versicolor
## 94 5.0 2.3 3.3 1.0 versicolor
## 95 5.6 2.7 4.2 1.3 versicolor
## 96 5.7 3.0 4.2 1.2 versicolor
## 97 5.7 2.9 4.2 1.3 versicolor
## 98 6.2 2.9 4.3 1.3 versicolor
## 99 5.1 2.5 3.0 1.1 versicolor
## 100 5.7 2.8 4.1 1.3 versicolor
## 101 6.3 3.3 6.0 2.5 virginica
## 102 5.8 2.7 5.1 1.9 virginica
## 103 7.1 3.0 5.9 2.1 virginica
## 104 6.3 2.9 5.6 1.8 virginica
## 105 6.5 3.0 5.8 2.2 virginica
## 106 7.6 3.0 6.6 2.1 virginica
## 107 4.9 2.5 4.5 1.7 virginica
## 108 7.3 2.9 6.3 1.8 virginica
## 109 6.7 2.5 5.8 1.8 virginica
## 110 7.2 3.6 6.1 2.5 virginica
## 111 6.5 3.2 5.1 2.0 virginica
## 112 6.4 2.7 5.3 1.9 virginica
## 113 6.8 3.0 5.5 2.1 virginica
## 114 5.7 2.5 5.0 2.0 virginica
## 115 5.8 2.8 5.1 2.4 virginica
## 116 6.4 3.2 5.3 2.3 virginica
## 117 6.5 3.0 5.5 1.8 virginica
## 118 7.7 3.8 6.7 2.2 virginica
## 119 7.7 2.6 6.9 2.3 virginica
## 120 6.0 2.2 5.0 1.5 virginica
## 121 6.9 3.2 5.7 2.3 virginica
## 122 5.6 2.8 4.9 2.0 virginica
## 123 7.7 2.8 6.7 2.0 virginica
## 124 6.3 2.7 4.9 1.8 virginica
## 125 6.7 3.3 5.7 2.1 virginica
## 126 7.2 3.2 6.0 1.8 virginica
## 127 6.2 2.8 4.8 1.8 virginica
## 128 6.1 3.0 4.9 1.8 virginica
## 129 6.4 2.8 5.6 2.1 virginica
## 130 7.2 3.0 5.8 1.6 virginica
## 131 7.4 2.8 6.1 1.9 virginica
## 132 7.9 3.8 6.4 2.0 virginica
## 133 6.4 2.8 5.6 2.2 virginica
## 134 6.3 2.8 5.1 1.5 virginica
## 135 6.1 2.6 5.6 1.4 virginica
## 136 7.7 3.0 6.1 2.3 virginica
## 137 6.3 3.4 5.6 2.4 virginica
## 138 6.4 3.1 5.5 1.8 virginica
## 139 6.0 3.0 4.8 1.8 virginica
## 140 6.9 3.1 5.4 2.1 virginica
## 141 6.7 3.1 5.6 2.4 virginica
## 142 6.9 3.1 5.1 2.3 virginica
## 143 5.8 2.7 5.1 1.9 virginica
## 144 6.8 3.2 5.9 2.3 virginica
## 145 6.7 3.3 5.7 2.5 virginica
## 146 6.7 3.0 5.2 2.3 virginica
## 147 6.3 2.5 5.0 1.9 virginica
## 148 6.5 3.0 5.2 2.0 virginica
## 149 6.2 3.4 5.4 2.3 virginica
## 150 5.9 3.0 5.1 1.8 virginica
## Predação
## 1 0.672
## 2 0.011
## 3 0.362
## 4 0.425
## 5 0.195
## 6 0.112
## 7 0.099
## 8 0.224
## 9 0.891
## 10 0.977
## 11 0.951
## 12 0.334
## 13 0.970
## 14 0.525
## 15 0.522
## 16 0.817
## 17 0.878
## 18 0.413
## 19 0.318
## 20 0.455
## 21 0.881
## 22 0.034
## 23 0.583
## 24 0.494
## 25 0.860
## 26 0.349
## 27 0.706
## 28 0.671
## 29 0.818
## 30 0.918
## 31 0.712
## 32 0.332
## 33 0.993
## 34 0.436
## 35 0.426
## 36 0.582
## 37 0.322
## 38 0.543
## 39 0.777
## 40 0.756
## 41 0.110
## 42 0.329
## 43 0.724
## 44 0.801
## 45 0.823
## 46 0.675
## 47 0.884
## 48 0.509
## 49 0.365
## 50 0.988
## 51 0.171
## 52 0.088
## 53 0.611
## 54 0.830
## 55 0.058
## 56 0.572
## 57 0.931
## 58 0.034
## 59 0.728
## 60 0.588
## 61 0.025
## 62 0.617
## 63 0.079
## 64 0.588
## 65 0.344
## 66 0.591
## 67 0.390
## 68 0.318
## 69 0.652
## 70 0.001
## 71 0.190
## 72 0.590
## 73 0.074
## 74 0.351
## 75 0.353
## 76 0.596
## 77 0.093
## 78 0.516
## 79 0.215
## 80 0.955
## 81 0.471
## 82 0.189
## 83 0.490
## 84 0.110
## 85 0.232
## 86 0.527
## 87 0.009
## 88 0.589
## 89 0.839
## 90 0.174
## 91 0.390
## 92 0.438
## 93 0.491
## 94 0.457
## 95 0.546
## 96 0.293
## 97 0.842
## 98 0.478
## 99 0.554
## 100 0.458
## 101 0.421
## 102 0.827
## 103 0.441
## 104 0.180
## 105 0.344
## 106 0.242
## 107 0.373
## 108 0.556
## 109 0.713
## 110 0.345
## 111 0.373
## 112 0.208
## 113 0.588
## 114 0.220
## 115 0.362
## 116 0.115
## 117 0.629
## 118 0.216
## 119 0.115
## 120 0.663
## 121 0.823
## 122 0.102
## 123 0.592
## 124 0.260
## 125 0.082
## 126 0.694
## 127 0.051
## 128 0.885
## 129 0.860
## 130 0.008
## 131 0.917
## 132 0.158
## 133 0.751
## 134 0.922
## 135 0.924
## 136 0.565
## 137 0.402
## 138 0.219
## 139 0.272
## 140 0.105
## 141 0.028
## 142 0.341
## 143 0.101
## 144 0.295
## 145 0.750
## 146 0.616
## 147 0.562
## 148 0.077
## 149 0.163
## 150 0.817O próximo passo é criar um modelo. Usaremos a função nativa
lm(), selecionando a variável Predação como
variável resposta e as demais (.) como preditoras. Em
seguidas, vamos aplicar a função nativa summary() para
extrair as estatísticas do modelo.
modelo <- lm(Predação ~ ., data = iris_pred %>% dplyr::select(dplyr::where(is.numeric)))
modelo %>%
summary()##
## Call:
## lm(formula = Predação ~ ., data = iris_pred %>% dplyr::select(dplyr::where(is.numeric)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52639 -0.23185 -0.00145 0.22223 0.51027
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.47101 0.25748 1.829 0.0694 .
## Sepal.Length -0.04245 0.07246 -0.586 0.5589
## Sepal.Width 0.07404 0.07503 0.987 0.3254
## Petal.Length 0.06307 0.07146 0.883 0.3789
## Petal.Width -0.18009 0.11873 -1.517 0.1315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2754 on 145 degrees of freedom
## Multiple R-squared: 0.07036, Adjusted R-squared: 0.04472
## F-statistic: 2.744 on 4 and 145 DF, p-value: 0.03077O modelo se apresentou como significativo (F(6,
143) = 2,188, p = 0,047), mas com um poder de explicação
baixo, explicando apenas 4% da variação (R2 = 0,04). Uma
coisa interessante desse modelo foi que nenhuma das variáveis preditoras
foram significativas (coeficientes estimados baixos, estatísticas
t baixas e p > 0,05). Esse efeito ocorre quando o modelo
sofre de multicolinearidade. Pelos índices de correlação previamente
feitos, sabemos que esse conjunto de dados apresentam
multicolinearidade. Vamos então testar se o VIF apresenta as mesmas
variáveis a serem removidas, utilizamos a função
car::vif(), do pacote car, utilizando nosso
modelo como os dados que entrarão na função. Vamos utilizar o ponto de
corte como 5 para o VIF.
library(car)## Carregando pacotes exigidos: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somemodelo %>% car::vif()## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 7.072722 2.100872 31.261498 16.090175Assim como nos índices de correlação, as variáveis que mais inflaram
o modelo foram as relacionadas às pétalas (VIFPetal.Length =
31,26; VIFPetal.Width = 16,09). O comprimento da sépala
(Sepal.Length) apresentou um valor de VIF maior que o ponto
de corte estabelecido (VIFSepal.Length = 7,07). Contudo, ela
está mais próxima do ponto de corte que as variáveis de comprimento e
largura da pétala. Então vamos remover apenas as variáveis sobre
pétalas, e ver se a multicolinearidade some por completo. Vamos criar
outro modelo, apenas com dados sobre sépalas.
modelo2 <- lm(Predação ~ Sepal.Length + Sepal.Width, data = iris_pred)
modelo2 %>%
summary()##
## Call:
## lm(formula = Predação ~ Sepal.Length + Sepal.Width, data = iris_pred)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53020 -0.21721 -0.03024 0.21384 0.56069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.59334 0.24091 2.463 0.0149 *
## Sepal.Length -0.06072 0.02754 -2.205 0.0290 *
## Sepal.Width 0.07582 0.05232 1.449 0.1494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2764 on 147 degrees of freedom
## Multiple R-squared: 0.05052, Adjusted R-squared: 0.0376
## F-statistic: 3.911 on 2 and 147 DF, p-value: 0.02214Note que agora o modelo está mais ajustado (F(2, 147) = 3,911, p = 0,022), apesar de explicar menos agora (R2 = 0,037). Contudo, note também que agora uma das variáveis, o comprimento da sépala, se tornou significativa para o modelo (t = -2,20, p = 0,029), apesar do coeficiente estimado baixo (-0,06 \(\pm\) 0,027). Vamos testar se a multicolinearidade sumiu com completo.
car::vif(modelo2)## Sepal.Length Sepal.Width
## 1.014016 1.014016Agora sim a multicolinearidade sumiu por completo. Ambas as variáveis não apenas estão dentro do ponto de corte estabelecido, como apresentam os mesmos valores (VIFSepal.Length, Sepal.Width = 1,01).
Um ponto importante sobre o VIF, é que o valores das variáveis do
modelo não mudam, não importa qual seja a variável resposta, desde que
sejam os mesmos valores das variáveis preditoras. Por causa disso, é
possível calcular o VIF sejam criar modelos, apenas utilizando suas
variáveis. Utilizando novamente uma função chamada
usdm::vif(), mas agora do pacote usdm, é
possível calcular o VIF apenas com os dados das variáveis.
library(usdm)## Carregando pacotes exigidos: terra## terra 1.7.71##
## Attaching package: 'terra'## The following object is masked from 'package:tidyr':
##
## extract##
## Attaching package: 'usdm'## The following object is masked from 'package:car':
##
## vifiris %>%
dplyr::select(dplyr::where(is.numeric)) %>%
usdm::vif()## Variables VIF
## 1 Sepal.Length 7.072722
## 2 Sepal.Width 2.100872
## 3 Petal.Length 31.261498
## 4 Petal.Width 16.090175Por mais que sim mais simples utilizar a função
usdm::vif(), ainda é recomendado utilizar a função
car::vif(), principalmente quando se quer avaliar os
efeitos da multicolinearidade nas estatísticas dos modelos.
Um cuidado muito importante na escolha entre índices de correlação e VIF é que o VIF gera apenas relações matemáticas de inflação de modelo, mas ele pode apontar que uma variável menos interessante, do ponto de vista biológico e ecológico, deveria ser removida, enquanto uma vaariável menos interessante deveria permanecer, enquanto que nos índices de correlação ainda é possível fazer esta análise de escolha. Por isso, o VIF é recomendo em duas principais situações, sendo a segunda a mais importante:
Conjunto muito grandes de variáveis, onde escolher, de forma manual, pode ser muito trabalhoso, e dificil de analisar multiplas combinações;
Principalmente em conjunto de variáveis que possuem a mesma importância biológica e ecológica para tentar explicar o fenômeno estudado, ou importâncias muito semelhantes, onde a remoção de qualquer variável não afetará o poder de explicação biológica e ecológica do modelo.